home *** CD-ROM | disk | FTP | other *** search
/ Meeting Pearls 1 / Meeting Pearls Vol 1 (1994).iso / installed_progs / gfx / lise2.1 / lise / src / expfit.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-03-31  |  6.7 KB  |  260 lines
  1. /*
  2.    Exponential Fit
  3. */
  4.  
  5. #include <stdio.h>
  6. #include <spec.h>
  7.  
  8. float *spc, *err, *tim, *usp1, *usp2, *uer1, *uer2;
  9. int   flg_v;
  10. int   bgrnd,xbeg,xend;
  11. float chisq,fit_a, fit_b;
  12.  
  13. help()
  14. {
  15.   printf("\n\nExponential Fit y = c + b * exp(a * x)\n");
  16.   printf("ExpFit file [-o file] [-ba n] [-xbeg n] [-xend n] [-last n]\n");
  17.   printf("options:\n");
  18.   printf("   file      spectrum to be fitted\n");
  19.   printf("   -o file   output spectrum instead of stdout\n");
  20.   printf("   -ba  n    background for spectrum (to avoid Background fit)\n");
  21.   printf("   -xbeg n   first channel to fit\n");
  22.   printf("   -xend n   last channel to fit\n");
  23.   printf("   -last n   takes arithmetic mean of last n channels as\n");
  24.   printf("             background. Default is, to find the background\n");
  25.   printf("             by fitting an exponential curve to the spectrum\n");
  26.   printf("   -v        verbose mode. Prints fitted parameters\n");
  27.   printf("   -vv       very verbose mode. Also Prints fitting progress\n");
  28.   printf("\n(C) Rainer Kowallik\n");
  29.   exit(0);
  30. }
  31.  
  32. /* ---------------------------------------------------------
  33.    Liner regression after logarithm of spectrum
  34.    This is only to calculate the chi**2, all
  35.    other values are lost.
  36.    UP TO NOW, THE ERRORBARS ARE NOT INCLUDED !
  37.    I HAVE FORGOTTEN MY OLD SUPERBASIC SOURCE AT HOME, SORRY.
  38.    --------------------------------------------------------- */
  39.  
  40. float linreg(spc,tim,err,n,a,z)
  41. float spc[],tim[],err[];
  42. int n,a,z;
  43. {
  44. int i,m,k;
  45. float ysoll,x,y,sumx,sumxx,sumy,sumyy,sumxy,nges,
  46.       steigung,offset,erg,
  47.       *ylog;
  48.  
  49.    ylog=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  50.  
  51.    nges = (float) (z-a);    /* Number of Elements */
  52.    sumx = 0.0;
  53.    sumxx = 0.0;
  54.    sumy = 0.0;
  55.    sumyy = 0.0;
  56.    sumxy = 0.0;
  57.  
  58.    for(i=a;i<z;i++) {
  59.       x = tim[i];
  60.       y = spc[i];
  61.       y = y - ((float) n);
  62.       if(y<0.1) y=1.0;                 /* trap overflow on log */
  63.       y = (float) log((double) y);     /* linearization of exp function */
  64.       ylog[i] = y;                     /* store for later use */
  65.       sumx = sumx + x;
  66.       sumxx = sumxx + (x*x);
  67.       sumy = sumy + y;
  68.       sumyy = sumyy + (y*y);
  69.       sumxy = sumxy + (x*y);
  70.    }
  71.    sumx = sumx / nges;
  72.    sumy = sumy / nges;
  73.    sumxx = sumxx - (nges * sumx * sumx);
  74.    sumxy = sumxy - (nges * sumx * sumy);
  75.    sumyy = sumyy - (nges * sumy * sumy);
  76.    steigung = sumxy / sumxx;
  77.    offset = sumy - (steigung * sumx);
  78.  
  79.    /* calculating chi**2 */
  80.  
  81.    erg=0.0;
  82.    for(i=a;i<z;i++) {
  83.       x = tim[i];
  84.       y = ylog[i];
  85.       ysoll = (steigung * x) + offset;
  86.       y = y - ysoll;
  87.       erg = erg + (y * y); 
  88.    }
  89.  
  90.    fit_a = steigung;
  91.    fit_b = offset;
  92.  
  93.    free(ylog);
  94.    return(erg);
  95. }
  96.  
  97. float sd(x,y)
  98. float x,y;
  99. {
  100. if(y!=0.0) return(x/y);
  101. return(0.0);
  102. }
  103.  
  104. main(argc,argv)
  105. int argc;
  106. char *argv[];
  107. {
  108. int  n,m,i,max,flg_ba,lastn;
  109. char z[80],comment[80],nam[80];
  110. float x,alpha,sum1,sum2;
  111. FILE *fp;
  112.  
  113.    flg_ba=0;            /* defaults to exponential fit */
  114.    bgrnd = -1;
  115.    flg_v=FALSE;
  116.    xbeg = 0;
  117.    xend = -1;
  118.  
  119.    if(argc < 2) {help(); exit(0);}
  120.    strcpy(nam,argv[1]); 
  121.    if(checkopt(argc,argv,"-ba",z)) {
  122.       bgrnd=atoi(z); flg_ba=1;
  123.    }
  124.    if(checkopt(argc,argv,"-xbeg",z)) {
  125.       xbeg = atoi(z);
  126.    }
  127.    if(checkopt(argc,argv,"-xend",z)) {
  128.       xend = atoi(z);
  129.    }
  130.    if(checkopt(argc,argv,"-last",z)) {
  131.       flg_ba=2; lastn = atoi(z);
  132.    }
  133.    if(checkopt(argc,argv,"-v",z)) flg_v=1;
  134.    if(checkopt(argc,argv,"-vv",z)) flg_v=2;
  135.  
  136.    spc=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  137.    err=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  138.    tim=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  139.    usp1=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  140.    usp2=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  141.    uer1=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  142.    uer2=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
  143.  
  144. /* -----------------------------------------------
  145.    reading spectrum 
  146.    ----------------------------------------------- */
  147.  
  148.    max=readspec(nam,usp1,uer1,tim,comment);
  149.    if(xend < 0) xend = max - 1;
  150.  
  151. /* -----------------------------------------------
  152.    generating error spectrum 
  153.    ----------------------------------------------- */
  154.    for(n=0;n<max;n++) {
  155.       uer1[n] = sqrt(usp1[n]);
  156.    }
  157.  
  158. /* --------------------------------------------
  159.    finding correct background
  160.    -------------------------------------------- */
  161.  
  162.    switch(flg_ba) {
  163.    case 0:
  164.       bgrnd = expbafit(usp1,tim,uer1,xbeg,xend);
  165.       chisq = linreg(usp1,tim,uer1,bgrnd,xbeg,xend);
  166.       break;
  167.    case 1: /* allready done when parsing arguments */
  168.       chisq = linreg(usp1,tim,uer1,bgrnd,xbeg,xend);
  169.       break;
  170.    case 2:
  171.       bgrnd = sumspc(usp1,max-lastn,max) / lastn;
  172.       chisq = linreg(usp1,tim,uer1,bgrnd,xbeg,xend);
  173.       break;
  174.    }
  175.  
  176.    fit_b = exp(fit_b);
  177.  
  178.    if(flg_v) {
  179.       printf("CHI**2 = %f\n",chisq);
  180.       printf("background = %d\n",bgrnd);
  181.       printf("b (t0)     = %f\n",fit_b);
  182.       printf("a (Tau)    = %f\n",-1.0/fit_a);
  183.    }
  184.  
  185. /* --------------------------------------------
  186.    calculating Fit curve
  187.    -------------------------------------------- */
  188.    
  189.    for(n = 0; n < max; n++) {
  190.       spc[n] = bgrnd + (fit_b * exp(fit_a * tim[n]));
  191.    }
  192.  
  193.  
  194. /* -------------------------------------------
  195.    now we can go and save the fit spectrum 
  196.    ------------------------------------------- */
  197.  
  198.    strcpy(z,comment);
  199.    n=instr("|",comment);
  200.    if(n>0) midstr(z,comment,0,n-1);
  201.    strcat(z," | ExpFit");
  202.    writespec("",spc,err,max,2,z);
  203.    free(uer1); free(uer2); free(usp1); free(usp2);
  204.    free(spc); free(err); free(tim);
  205.    exit(0);
  206. }
  207.  
  208. sumspc(spc,a,z)
  209. float spc[];
  210. int a,z;
  211. {
  212. int i,erg;
  213.  
  214.    erg=0;
  215.    for(i=a;i<z;i++) erg=erg+spc[i];
  216.    return(erg);
  217. }
  218.  
  219.    
  220. /* ---------------------------------------------------------
  221.    Fitting background to exponential function.
  222.    This is done, by calculating the chi**2 with
  223.    linear regression for the logarithm of the spectrum.
  224.    We than look for a minimum of chi**2 for different
  225.    backgrounds.
  226.    UP TO NOW, THE ERRORBARS ARE NOT INCLUDED !
  227.    I HAVE FORGOTTEN MY OLD SUPERBASIC SOURCE AT HOME, SORRY.
  228.    --------------------------------------------------------- */
  229. expbafit(spc,tim,err,a,z)
  230. float spc[],tim[],err[];
  231. int a,z;
  232. {
  233. int i,n,erg,bamin,bamax,babest,inc,xraster,xdepth;
  234. float chibest;
  235.  
  236. xraster = 10;
  237. xdepth = 5;
  238. chibest = 1E10;
  239. bamax = sumspc(spc,z-10,z) / 10 ; /* get maximum background */
  240. bamin = bamax / 2;                /* get minimum background */
  241. for(i=1;i<xdepth;i++) {
  242.    inc = (bamax - bamin) / xraster;
  243.    if(inc<1) break;
  244.    n = bamin; while(n<=bamax) {
  245.       chisq = linreg(spc,tim,err,n,a,z); /* just calculate chi**2 */
  246.       if(chisq < chibest) {
  247.      chibest = chisq;
  248.      babest = n;
  249.      if(flg_v > 1) {
  250.      printf("new best: background = %d   chi**2 = %f\n",babest,chibest);
  251.      }
  252.       }
  253.       n = n + inc;
  254.    }
  255.    bamin = babest - inc;
  256.    bamax = babest + inc;
  257. }
  258. return(babest);
  259. }
  260.